…plementation

…lcuation in DCON)

…PEC-DCON

- Consolidated main program into a single `Main` function for better readability and modularity.
- Introduced `PentrcControl`, `PentrcInternal`, and `PentrcOutput` structures to encapsulate control parameters, internal state, and output configurations respectively.
- Created a new `Output.jl` file to manage output operations, including writing torque, orbit, and energy data to ASCII files.
- Improved error handling and output directory management.
- Updated documentation and comments for clarity and consistency.

- Updated module descriptions and added detailed function documentation for energy integration and torque calculations.
- Refactored special functions to include elliptic integrals and double factorial computations.
- Improved clarity and organization of code for better maintainability.

…d functions in Splines module

…management; implement spline_roots and spline_fit functions in CubicSpline

…RC into GeneralizedPerturbedEquilibrium

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

…States; add documentation page

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

…pline_roots with Roots.jl

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

…fix scoping, resolve duplicates

- energy.jl: Replace TaskLocalValue with plain Dict (no TaskLocalValues.jl dep)
- energy_integration.jl: Remove duplicate xintgrl_lsode stub (real impl in energy.jl),
  remove duplicate global Ref() definitions, add verbose keyword to tintgrl_grid
- pentrc.jl: Add includes for energy_integration.jl and special.jl
- Input.jl: Move `using DelimitedFiles` to top level, remove `using HDF5` from
  inside functions (already in pentrc.jl), remove shadowed local const declarations
- Utils.jl: Remove dead Utilities sub-module with broken params.jl include

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

…etic energy matrix producer

In `compute_kinetic_matrices_at_psi!`, the producer assigned `kwmat = complex(0, imag)`
and `ktmat = complex(real, 0)`, but the downstream FKG reduction in
`ForceFreeStates/Kinetic.jl` follows the Fortran convention where `kwmat` carries
pure-real values (from the "wmm" integrand, `rex=0, imx=1`) and `ktmat` carries
pure-imag values (from the "tmm" integrand, `rex=1, imx=0`), both emerging from the
-i/(2n) normalization. The swap corrupted the asymmetric kinetic adjoints
(`caat = cmat_kin - 2·ktmat[3]`, `bkmat = kwmat[2] + ktmat[2] + i·(...)`), producing
sign-flipped F̄, K̄, Ḡ in the kinetic ODE RHS.

Validated on Solovev kinetic regression (Re(et[1]) sign flips -115.7 -> +99.75)
and the DIIID kinetic example (mpsi=128, pow1 grid, hamada coords, delta_m=±8,
set_psilim_via_dmlim): Re(et[1]) = +35, 4 singular surfaces q=2..5, no pathology.
Ideal Solovev regression invariant to ~1e-13.

Also rename for clarity:
- calculate_gar_matrices! -> kinetic_energy_matrices_for_euler_lagrange!
- out_wmats/full_wmats    -> cmplx_kinetic_energy_mats

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

…int before FFT

Fortran fspline_fit_2 (math/fspline.f:293) uses f = fst%fs(:, 0:my-1, iq),
explicitly excluding the duplicated endpoint. Julia was FFT'ing all length(ys)
samples including the θ=0 ≡ θ=2π duplicate, biasing the DC coefficient by
~(f(0) − mean)/N. Small effect (~1% for the a10 kinetic Bk matrix at
ψ=0.30928), but principled — matches Fortran convention.

Detects the duplicate automatically via ys[end] − ys[1] ≈ 2π; if ys is already
non-duplicated, the FFT runs on all samples as before.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

…gral

The phase factor pl[i] = exp(-2πi·lnq·cum_wb_arr[i]/pl_denom) used a
left-Riemann cumulative sum of the bounce integrand, whereas Fortran
torque.F90:736-750 uses spline_fit + spline_int, which for smooth functions
produces a cubic-spline cumulative integral ≈ trapezoidal:

    Fortran bspl%fsi(j)/Δx ≈ g_1 + g_2 + ... + g_{j-1} + g_j/2
    Julia   cum_wb_arr[i]   = g_1 + g_2 + ... + g_{i-1}    (too large by g_{i-1}/2)

The discrepancy was a per-sample phase offset of ~π·lnq·g_i/pl_denom, small
per-point but concentrating at the m=0 column of the kinetic matrices
because m=0 has no oscillating Fourier basis factor to average the
sampling bias out.

Effect at ψ=0.30928, a10 kinetic example (Bk matrix):
- Re(Bk[:, m=0]) Julia/Fortran ratio: 10.5 → 1.00
- Full Bk rel_frob:                    4.30% → 1.60%
- Hk rel_frob:                         30-100% → 1.61%
- Ck, Dk, Ek rel_frob:                 2-5% → 0.5-1.9%

Fix: cum_wb_arr[i] = cum_wb − wb_integrand/2 (same for cum_wd_arr).

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

… integrals

Convert the three right-Riemann sums in _bounce_integrate (bj_integral and
the wmu_ba/wen_ba bounce-average loops) to explicit trapezoidal with 0.5
endpoint weights. The sums were formally right-Riemann over 2:ntheta,
bit-equivalent to trapezoidal only because the payload arrays (jvtheta,
wmu_mt, wen_mt) are zero at i=1 and i=ntheta from the 2:ntheta-1 population
loop. Writing the weights explicitly makes the quadrature self-correct if
the boundary handling ever changes.

Also remove cum_wd_arr: allocated and written since introduction
(163aef2) but never read. Mirrors a commented-out Fortran branch
(torque.F90:747-751) for magnetic-precession-dominated regime that neither
code activates.

Numerical effect: all 6 kinetic matrix rel_frob values at a10 ψ=0.30928
unchanged (Ak 2.82%, Bk 1.60%, Ck 1.87%, Dk 0.54%, Ek 0.77%, Hk 1.62%).
Tiny ~1e-6 residuals in kwmat Re and ktmat Im (noise from fwmm/ftmm
resonance-operator symmetry) go to exact zero.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

…t/GPEC into feature/pe-kf-interface

…ic kwmat/ktmat

Adds `integrate_pitch_gar_quadgk_wt` and `_pitch_gar_kernel_quadgk_wt!` that emit
BOTH the fwmm half (rex=0, imx=1 → Fortran kwmat) and the ftmm half (rex=1, imx=0
→ Fortran ktmat) in a single pitch integration, sharing one energy integration per
(λ, E). Returns a length-`2*nqty` packed buffer `[wmm | tmm]`.

The prior single-integration approach (rex=imx=1, then split kwmat=real, ktmat=imag)
gives correct forward sums kwmat+ktmat = full but WRONG adjoint combinations
kwmat-ktmat for non-Hermitian B_k, C_k, E_k blocks — because Fortran's kwmat and
ktmat are each genuinely complex (matching Fortran pentrc/torque.F90:842-847 rex/imx
assignments and pentrc/pitch.f90:376 integrand decomposition), not pure real / pure
imaginary. Per-surface matrix dumps against Fortran dcon/fourfit.F:1080-1082
(`kwmat_l`, `ktmat_l`) confirm element-by-element agreement with the dual-output
convention.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

…-triangle reconstruction

Rewires `kinetic_energy_matrices_for_euler_lagrange!` and
`compute_kinetic_matrices_at_psi!` to populate `kwmat` and `ktmat` directly via the
new `integrate_pitch_gar_quadgk_wt` dual-output kernel, replacing the prior post-hoc
`complex(real(full), 0)` / `complex(0, imag(full))` split.

For A, D, H Hermitian-outer-product blocks (wz†wz, wx†wx, wy†wy) stored as upper
triangles, the lower-triangle reconstruction now uses different mirrors per half:

  kwmat[j,i] =  conj(kwmat[i,j])   (Hermitian)
  ktmat[j,i] = -conj(ktmat[i,j])   (anti-Hermitian)

Derivation: S_w = complex(0, imag(xint)) is pure imaginary so conj(S_w) = -S_w,
whereas S_t = complex(real(xint), 0) is pure real so conj(S_t) = S_t. Combined with
conj(factor) = -factor (factor = -i/(2n)), these mirrors recover Fortran's
independently-computed (j,i) slots exactly. Implemented via a `mirror_sign`
parameter to a shared `_assemble_hermitian!` closure.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

…itian kinetic A

The previous implementation computed `aamat = (amat_lu \ amat_kin)'` which equals
`(A_kin⁻¹ · A_kin)' = I` exactly. This silently zeroed `umat_diff = I - aamat` and
DROPPED the `im·psio_over_n · umat_diff · b1mat` term from `paat`, the
`im·psio_over_n · aamat · bkmat` term from `r1mat`, and the
`im·psio_over_n · umat_diff · cmat_kin` term from `r2mat` whenever `amat_kin` was
non-Hermitian.

Fortran dcon/sing.f:1004-1008 computes

    temp2 = amat_kin
    zgbtrs("C", ..., amatlu, temp2)        ! temp2 = A_kin⁻ᴴ · amat_kin
    aamat = CONJG(TRANSPOSE(temp2))        ! aamat = amat_kin^H · A_kin⁻¹

which is NOT the identity when `amat_kin` is non-Hermitian. This was hidden while
the kinetic producer gave pure-real `kwmat` and pure-imag `ktmat` (amat_kin stayed
Hermitian); exposing it required the correct dual-output kernel (prior commits) that
makes amat_kin genuinely non-Hermitian via the anti-Hermitian ktmat contribution.

Fix: `aamat = (amat_lu' \ amat_kin)'` — the `amat_lu'` backslash path solves Aᴴ·x = b,
giving `temp = A_kin⁻ᴴ · amat_kin` which after adjoint becomes the Fortran form.

With all three kinetic fixes in place (dual-output kernel, anti-Hermitian ktmat
mirror, and this aamat correction), the DIIID_kinetic_example DCON benchmark
reaches Fortran agreement:

  Re(et[1]): Julia 1.0481, Fortran 1.0507, err 0.25%
  Im(et[1]): Julia -0.3012, Fortran -0.3013, err 0.01%

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

…ts_prim

Clarifies that Julia stores the primitive (pre-Schur-reduction) geometric forms
D = χ₁·(g23 + q·g33·m/n) and the analogous E. The `_prim` suffix follows the
existing `fmats_prim` convention.

Previously, a reader diffing against Fortran fourfit would find that Fortran
saves two distinct splines — `dbats/ebats` (primitive) and `dmats/emats`
(kinetic-block overwrites with χ₁²·(g22+2q·g23+q²·g33) and related) — whereas
Julia's field named `dmats` actually holds the primitive form. The overwritten
form is only consumed by Fortran's alternate on-demand singular-surface Schur
path (`fkg_kmats_flag=.false.` in sing.f), which Julia does not implement.
Renaming makes the intent explicit and avoids a silent naming trap.

Pure rename; no numerical change.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

…rface

# Conflicts:
#	src/Equilibrium/DirectEquilibrium.jl

…lytic resonance poles

Port the QuadGK energy-integration robustness from PR #220 into the
KineticForces energy integrator, with a corrected pole decomposition.

- Map x = E/T to u = 1-exp(-x) on [0,1), absorbing the Maxwellian weight
  into du, so the integral covers [0,infinity) without an xmax cutoff.
- Remove each resonance pole analytically by a Sokhotski-Plemelj
  decomposition: subtract R/(u-u_pole) and add back the elementary
  integral R*(log(1-u_pole)-log(-u_pole)). The same complex pole u_pole
  (x_pole = x_res - i*nu/Omega') is used in both the subtraction and the
  add-back, so the decomposition is exact and converges for nu > 0 as well
  as nu -> 0 (PR #220 subtracted at the real u_res but added at the
  complex u_pole, which diverges for collisional cases).
- Collisionless limit uses the causal Sokhotski-Plemelj branch.
- CGL has no resonance denominator: integrate the physical x-space
  integrand directly over [0,infinity) via QuadGK.
- ximag is accepted for signature compatibility but is now unused.

Adds find_resonance_energies and unit tests: resonance-root cases, a
collisional cross-check against direct x-space quadrature, and a
collisionless-limit continuity check. 110 kinetic unit tests pass.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

…rface

…e resonances

When kinetic_source="calculated" runs the full KF torque pipeline, the pitch
loop can hand the energy integrator parameter sets whose resonance polynomial
has a root at very large or very small x. Two corresponding NaN traps appear
in the Sokhotski-Plemelj analytic pole contribution:

- x_res > 700: exp(-x_res) underflows to zero, R becomes 0+0im, and
  R*(log(1-u_pole) - log(-u_pole)) evaluates as 0*(-Inf+iπ) = NaN. The
  resonance's true contribution is identically zero (the Maxwellian weight
  is below machine precision) — skip the resonance.
- nu_res / omega_prime = Inf (harmonic ν overflows at tiny x_res, or
  omega_prime ≈ 1e-30): exp(-x_pole) at infinite imaginary argument returns
  NaN. The resonance is infinitely collisionally broadened — no localized
  pole, no extraction needed; skip and let QuadGK integrate the smooth
  integrand directly.

Verified by the regression-harness solovev_kinetic_calculated case (which
exercises the calculated-source KF pipeline that the "fixed" fullruns do
not): pre-port et[1] = 19.914 - 0.628i, post-fix et[1] = 19.908 - 0.658i
(0.03% / 5% relative). 110 kinetic unit tests still pass.

Also updated:
- test/test_data/regression_solovev_kinetic_calculated/gpec.toml: replace
  removed dmlim/set_psilim_via_dmlim keys with develop's psiedge edge-scan
  band so the calculated-source case is runnable.
- test/runtests_fullruns.jl: refresh ex4 (Solovev kinetic multi-n) baseline
  to the develop-consistent eigenvalue 0.22325 — the prior -0.01248 dates
  from before develop's edge-scan and periodic-theta-endpoint evolution.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

…schema

The DIIID kinetic benchmark's auto-generated gpec.toml used the deprecated
set_psilim_via_dmlim and dmlim keys, which develop removed (replaced with
the psiedge edge-scan band). The new FFS rejects unknown keys, so the
benchmark errored out before any compute. Drop the deprecated pair; the
existing psiedge = 1.00 (no edge-scan truncation) is the closest current
equivalent.

A/B verified with the new u-substitution energy integrator on the
merged tree: same DIIID benchmark numbers as the pre-port integrator
(FGAR T = 1.805 N·m vs Fortran 1.780, 1.39% err; FGAR dW = 0.160 vs
0.137, 16.94% err). The dW drift from PR #224's body (0.141 → 0.160)
is caused by the develop merge upstream of KF (most likely the
periodic-theta-endpoint fix changing bounce-data evaluation), not by
this port — the energy integrator is faithful.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

KineticForces — Complete PE→KF pipeline and validate DIIID-PENTRC benchmark